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Abstract. We present hydrodynamical calculations of radiative shocks with low Mach numbers and find that the 
well-known global overstability can occur if the temperature exponent (a) of the cooling is sufficiently negative. 
We find that the stability of radiative shocks increases with decreasing Mach number, with the result that M = 2 
shocks require a < —1.2 in order to be overstable. Such values occur within a limited temperature range of many 
cooling curves. We observe that Mach numbers of order 100 are needed before the strong shock limit of a cr ~ 0.4 
is reached, and we discover that the frequency of oscillation of the fundamental mode also has a strong Mach 
number dependence. We find that feedback between the cooling region and the cold dense layer (CDL) further 
downstream is a function of Mach number, with stronger feedback and oscillation of the boundary between the 
CDL and the cooling region occuring at lower Mach numbers. This feedback can be quantified in terms of the 
reflection coefficient of sound waves, and in those cases where the cooling layer completely disappears at the end 
of each oscillation cycle, the initial velocity of the waves driven into the upstream pre-shock flow and into the 
downstream CDL, and the velocity of the the boundary between the CDL and the cooling layer, can be understood 
in terms of the solution to the Riemann problem. An interesting finding is that the stability properties of low 
Mach number shocks can be dramatically altered if the shocked gas is able to cool to temperatures less than the 
pre-shock value (i.e. when \ < 1, where \ is the ratio of the temperature of the cold dense layer to the pre-shock 
temperature). In such circumstances, low Mach number shocks have values of a cr which are comparable to values 
obtained for higher Mach number shocks when x = 1- F° r instance, a CI = —0.1 when M = 2 and \ = 0.1, 
comparable to that when M = 10 and \ = 1. Thus, it is probable that low Mach number astrophysical shocks 
will be overstable in a variety of situations. We also explore the effect of different assumptions for the initial 
hydrodynamic set up and the type of boundary condition imposed downstream, and find that the properties 
of low Mach number shocks are relatively insensitive to these issues. The results of this work are relevant to 
astrophysical shocks with low Mach numbers, such as supernova remnants (SNRs) immersed in a hot interstellar 
medium (e.g., within a starburst region), and shocks in molecular clouds, where time-dependent chemistry can 
lead to overstability. 
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1. Introduction 

Radiative shocks were first shown to be susceptible to 
a cooling instability by Falle (1975) and McCray et 
al. ©75), and to oscillations from a global overstability 
by Langer et al. (1981), and have since been extensively 
examined (see Pittard et al. 2003 and references therein) . 
A central finding is that the stability of a radiative shock 
deteriorates as the Mach number increases (Strickland & 
Blondin 1995), since the higher density of the cold gas 
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layer behind radiative shocks of higher Mach number acts 
increasingly like a reflecting wall, and reflects incident 
waves with virtually the same amplitude. In contrast, in- 
cident waves are reflected with much decreased amplitude 
in lower Mach number shocks since the lower density cold 
gas layer acts more like a cushion. 

In most investigations of the radiative overstability, 
the temperature dependent cooling function is approxi- 
mated as A(T) = AoT", with the rate of energy loss per 
unit volume, E — n 2 A(T). The radiative overstability is 
known to depend upon the value of a: if a is greater than 
some critical value, a cr , the system is stable to the fun- 
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damental mode of oscillation (though it may still be un- 
stable to higher overtones). Previous numerical work has 
shown that for high Mach numbers shocks a cr w 0.4 (e.g., 
Imamura et al. 119841 Strickland & Blondin 1995), in good 
agreement with the linear stability analysis of Chevalier 
& Imamura |T2H2|- The first overtone mode is stabilized 
when a > 0.8. Further work has revealed that non-radial 
transverse modes may be unstable at values of a which are 
stable to radial modes (cf. Bertschinger I198T>)) . and that 
non-radial effects are greater at higher densities, perhaps 
because of additional types of instabilities (see Blondin et 

To date, most investigations of the global overstability 
of radiative shocks have been concerned with high Mach 
number shocks (M > 5). While simulations of decelerating 
radiative shocks in SNRs (e.g., Kimoto & Chernoff 119971 
Blondin et al. H998[) and wind-blown bubbles (Walder & 
Folini I1996f) "scan" a wide range of Mach numbers, the 
ISM densities and temperatures considered in these works 
resulted in SNRs with slow (due to ISM density) low Mach 
number (due to ISM temperature) forward shocks, with 
post-shock temperatures (slightly higher than the ISM 
temperature) which lie in a range on the cooling curve 
where a is much larger than zero. Thus there is no insta- 
bility. However, instability may occur if the temperature 
of the ISM is larger, such that the cooling function then 
has a significantly steep temperature dependence, though 
there is no information in the current literature on the 
value of a cr for low Mach number shocks. The ISM tem- 
perature 1 assumed in existing papers is typically ~ 10 4 K. 

In this paper we therefore examine the overstability 
of low Mach shocks, extending previous numerical work 
to Mach numbers below 5. Since the boundary conditions 
imposed in numerical models of radiative shocks can also 
affect their stability (Strickland & Blondin 1995, Saxton 
2002}, we investigate the nature of the instability as the 
boundary conditions and initial conditions are varied. In 
Sec. |21 we summarize the properties of low Mach number 
radiative shocks. Hydrodynamical models of the oversta- 
bility are presented in Sec. |3] where we also examine how 
feedback between the cold, dense layer (CDL) downstream 
of the cooling layer, and the cooling region depends on the 
Mach number. Additional results obtained when we allow 
the shocked gas to cool to a temperature lower than that 
of the pre-shock flow are noted. We discuss applications 
of our work in Sec. 01 and Sec. summarizes the paper. 

2. Properties of low Mach number shocks 

The approximation that A(T) = AoT" is particularly ap- 
propriate for the low Mach number flows which are the 
subject of this work because the shocked gas cools through 

1 Usual practice is to set the floor-temperature of the simu- 
lation and Tism to the same value. Most of the simulations in 
this paper adopt this assumption, though the effect of letting 
the gas cool below the pre-shock temperature is examined in 
Sec. 13.41 where we find that there are dramatic changes in the 
resulting properties when the Mach number is low. 



Table 1. Properties of a radiative shock for gas with a 
ratio of specific heats, 7 = 5/3. 



M 



Ps/po pc/po T s /T 



1.4 


1.58 


3.27 


1.39 


2.0 


2.29 


6.67 


2.08 


3.0 


3.00 


15.0 


3.66 


5.0 


3.57 


41.7 


8.68 


10.0 


3.88 


167 


32.1 


20.0 


3.97 


667 


126 


40.0 


3.99 


2667 


501 



a relatively narrow range of temperatures. Some proper- 
ties of radiative shocks are tabulated in TablefJ] The ratio 
of post- to pre-shock density is given by 

Ps _ (7 + l)Af 2 
po (7 - 1)M 2 + 2 

the ratio of post- to pre-shock temperature is given by 
7, = [2 7 M 2 - (7 - 1)][( 7 - 1)M 2 + 2] 
To 



(1) 



(2) 



[(7 + 1)M]2 
and the ratio of final to pre-shock density is 

^ = 7M 2 (3) 
Po 

if the gas cools to its initial temperature, To. In Eqs. ^ 
|U the subscript "s" denotes quantities immediately be- 
hind the shock front, while the subscript "0" indicates 
pre-shock quantities. The subscript "c" indicates quan- 
tities associated with the CDL. One consequence of the 
reduced ratio of T s /Tq in low Mach number shocks is that 
a can obtain fairly extreme negative values when the lo- 
cal slope of the cooling curve becomes very steep. This is 
demonstrated in Figs. Q and [3 where we show two exam- 
ple cooling curves and the slope a. In both curves there 
are two distinct temperature ranges where a < —2, and in 
the cooling curve where collisional ionization equilibrium 
is enforced we find that a k, —4 when T « 3 x 10 5 K. We 
further note that time-dependent chemistry behind shocks 
in molecular clouds can also result in effective values of a 
which are extremely negative (Smith & Rosen 2003). 

The cooling length of a radiative shock can be param- 
eterized as 

w s mfcT s 



L c = x h L c 



Xh- 



(4) 



MoT"' 

where p, v and T are the fluid mass density, velocity, and 
temperature, respectively, rh is the mean mass per par- 
ticle (p = fan where n is the number density) and k is 
Boltzmann's constant. L c is defined as the distance be- 
hind the shock that a parcel of gas travels before cooling 
to the pre-shock temperature, at which point cooling is 
usually assumed to turn off. x\, is a factor of order unity. 
The structure of the shock and the value of x^ can be 
determined by the method noted in Strickland & Blondin 
(1995)). Values for L c and Xb are given in Tabled for the 
case where m, k and Ao are equal to unity. This assump- 
tion is also adopted for the numerical work which follows. 
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Fig. 1. Cooling curves for material of solar abundance as- 
suming collisional ionization equilibrium (solid) or non- 
equilibrium ionization (dashed). The former was calcu- 
lated using the MEKAL thermal plasma code (Mewe et 
al.EHEHI KaastraHHHSl distributed in XSPEC (vll.2.0), 
while the latter was taken from data in Sutherland & 
Dopita ifT^Mjl . These curves are normalized so that the 
net cooling rate per unit volume, E = n c n,iA(T), where 
n e is the electron number density and ni is the total num- 
ber density of all of the ions. 




-3 1 — ' ^ — 1 1 1 1 11 — ' ^ — ' 

10 4 10 5 10 6 10 7 10 1 
T(K) 



Fig. 2. The slope, cv, of the respective cooling curves in 
Fig.EKCIE - solid; NEI - dashed). 

3. Hydrodynamical Simulations 

In this section we describe the results obtained with the 
hydrodynamical code VH-1 (see Blondin et al. 1990). VH- 
1 uses a Lagrangian formulation of the piecewise parabolic 
method (PPM; see Colella & Woodward ITMljl . where the 
spatial representation of the fluid variables is third order 
accurate. After each timestep the conserved quantities are 
remapped onto the original grid. Recent work in which 
this scheme was used for investigations of the overstability 
of radiative shocks includes Strickland & Blondin (1995), 
Sutherland et al. QMSty . and Pittard et al. ifl^fMl l27)M|) . 

Radiative cooling is implemented via operator split- 
ting and an unconditionally stable implicit scheme (see 
Strickland & Blondin 1995 ). At unresolved phases between 
hot gas and denser cold gas, the rate of energy loss of a 



Table 2. Cooling length of a radiative shock for gas with 
a ratio of specific heats, 7 = 5/3, and cooling exponent, 
a = 0. L' c in appropriate units can be obtained by mul- 
tiplying by the numerical values of m&/Ao (note that in 
this work we also specify that p — 1 and vq = 1). There 
is convergence in the value of x\, as M increases. 



M 


L' e 


x h 


L c 


1.40 


0.171 


0.545 


9.294E-02 


2.00 


5.966E-02 


0.685 


4.087E-02 


3.00 


2.716E-02 


0.748 


2.032E-02 


5.00 


1.633E-02 


0.769 


1.256E-02 


10.0 


1.278E-02 


0.776 


9.917E-03 


20.0 


1.198E-02 


0.777 


9.307E-03 


40.0 


1.178E-02 


0.778 


9.167E-03 



given hydrodynamical cell due to radiative cooling is re- 
stricted to the minimum from the neighbouring cells. This 
procedure allows the rapid cooling layer at the rear of 
radiative shocks to be "resolved" with relatively few hy- 
drodynamical cells, and provides the advantage of a sig- 
nificant speed-up to computational times. Resolution tests 
without this restriction show that the oscillation frequency 
and amplitude of the overstable shock converge towards 
the results obtained with this restriction as the numerical 
resolution is increased (see also Sutherland et al. 2003 ) . 
We have performed additional tests to check that the be- 
haviour of the system converges with increasing numerical 
resolution when a — a cr . 

While a substantial body of work on the overstability 
of high Mach number shocks exists, it is difficult to make 
comparisons between results since the adopted bound- 
ary conditions often differ. It has been shown that the 
boundary conditions can play an important role in de- 
termining the stability of the shock (see Strickland & 
Blondin 1995). Therefore, we have examined the effects 
of 3 different types of initial and boundary conditions. 
Results from each are presented in turn, where we prevent 
shocked gas from cooling below the pre-shock tempera- 
ture. Simulations where we allow the gas to cool further 
are presented in Sec. 13.41 

3.1. Impulsive shock generation 

In this set of simulations we run a supersonic flow into a 
pre-existing cold, dense layer of g shown schemati- 

cally in Fig. The density and velocity of the supersonic 
flow are set to unity, and its pressure and temperature are 
set to I/7M 2 . The density of the cold layer is set to 7M 2 , 
and its thermal pressure is set to unity (i.e. equal to the 
ram pressure of the oncoming flow). The velocity of gas 
within the cold layer is set to I/7A/ 2 . The total pressure 
(thermal plus ram) within the supersonic flow is equal to 
that within the cold dense slab. An inflow boundary ex- 
ists in the supersonic flow, while an outflow boundary is 
specified at the edge of the CDL. All simulations are one- 
dimensional and use 7 = 5/3. 
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Fig. 3. Schematic of the grid set up for implusive shock 
generation. Supersonic flow with density po and flow speed 
Vq collides with a subsonic CDL with density jM 2 po and 
flow speed vq/^M 2 . The upstream boundary condition is 
inflow, while the downstream boundary condition is out- 
flow. All simulations are 1-D. 



Table |21 also reveals that a M = 5 radiative shock has 
a value of a CT which is still far from that which occurs in 
the strong shock limit (a cr « 0.4, cf. Chevalier & Imamura 
1982). Calculations reveal that this limit is not reached 
until Mach numbers of ss 100, with shocks of M = 10, 20, 
and 40 having a CI ~ —0.1, 0.2 and 0.3, respectively. 
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We set up the grid such that the width of each hydro- 
dynamical cell is constant in the supersonic flow and for 
several cooling lengths (L c ) into the CDL, after which the 
width of each cell increases by 3% relative to its neighbour. 
This enables the outflow boundary to be positioned so far 
downstream that it does not affect the evolution of the 
overstability (i.e. the shock driven into the CDL does not 
reach the outflow boundary by the time we stop our sim- 
ulations). Our simulations typically contain ~ 300 — 500 
cells within the cooling length of the shock. 

Upon starting the simulation, a radiative cooling layer 
is impulsively generated as the supersonic flow plows into 
the dense cool layer. This set up is similar to that adopted 
by Sutherland et al. l|2003fl but with the difference that 
in our work the outflow boundary is much further down- 
stream. In Fig. 0| we show the time-dependent behaviour 
of simulations with M = 1.4 and a = —2.5, —2.0, and 
— 1.5. It is clear from Fig. 0] that the overstability be- 
comes damped as a increases, but exists if a is small 
enough. When M — 1.4, we find damping is achieved with 
a > —1.8. Comparison with Fig. reveals that the value 
of a required for the overstability of low Mach number 
flows occurs within certain limited temperature ranges. 
In Fig. [3] we show the time-dependent behaviour of sim- 
ulations with a = —1.5 and M — 1.4, 2, 3, and 5. It is 
clear that radiative shocks become increasingly suscepti- 
ble to overstability as M increases, a finding which is in 
agreement with earlier work. In Table|3|we list the critical 
value of a for damping of the fundamental mode. Note 
that nonlinear effects may contribute to these values (see, 
e.g., Strickland & Blondin fl995|l . We discuss this further 
in connection with a linear stability analysis of low Mach 
number shocks (Pittard et al. in preparation) . 




Fig. 4. Time-space diagrams of the density evolution of 
a 1-dimensional Mach 1.4 radiative shock with different 
cooling exponents a. Supersonic flow enters the grid from 
the top, and the cooled postshock gas flows off the grid 
at the bottom. Lighter shades indicate lower densities. 
Distances are marked in units of L c (the value of this 
is different in each panel - see Eq. 0] and Table |2J) while 
time is shown in units of L c /vq (vq = 1). The value of a 
is —2.5, —2.0, and —1.5 in the top, middle, and bottom 
panels respectively. 

Power spectra of the shock position of the simulations 
shown in Fig. 0] reveal that the dimensionless angular fre- 
quency of the fundamental mode, u>{, is relatively insensi- 
tive to a, being 0.665, 0.690, and 0.734 for a of -2.5, -2.0 
and —1.5 respectively 2 . This trend is consistent with ear- 
lier simulations of higher Mach number shocks (Strickland 
& Blondin lTj|9 5| ■ However, larger changes in ojf occur 
when a is kept constant and M is varied. This is illus- 
trated in Fig. [fj] where power spectra of the simulations 



2 For ease of comparison with previous work, uji is defined 
in units of vo/xb , where vo is dimensionless and equal to unity 
and where x b is also dimensionless. 
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Fig. 5. Time-space diagrams of the density evolution of 
a 1-dimensional radiative shock with cooling exponents 
a = —1.5 and different Mach numbers. M is varied from 
1.4, 2, 3 and 5 from the top to bottom panels. Note that 
the axis and density scaling changes from panel to panel, 
unlike in Fig. U where the scaling is kept constant. 

shown in Fig.|5]are displayed. Table[3]also notes the value 
of oj[ as a function of M for shocks with a = —1.5. As M 
increases, Wf decreases. 

3.2. Alternative set ups 

One alternative to the work in the previous section is to 
impose, at the beginning of the simulation, a steady-state 
shock structure in the middle of the grid. The oversta- 
bility may then be excited from weak perturbations pro- 
duced by the numerical mapping of the steady solution 
onto the computational grid. While the initial conditions 
are clearly different, the boundary conditions are identical 
to the previous case (that is inflowing pre-shock gas at the 
upstream boundary, and outflowing, cold, compressed gas 



Table 3. Critical value of a for damping of the funda- 
mental mode, and angular frequency (when a = —1.5) 
as a function of the Mach number. For a > cv cr the fun- 
damental mode is damped, though the shock may still 
be unstable to higher order modes, such as the first over- 
tone (cf. Chevalier & Imamura 1982, Strickland & Blondin 

mm . 



M 




LUf (a = —1.5) 


1.40 


-1.8 


0.73 


2.00 


-1.2 


0.51 


3.00 


-0.7 


0.32 


5.00 


-0.4 


0.22 



at the downstream boundary) . Calculations based on this 
method are presented in Strickland & Blondin (1995) and 
Pittard et al. QXKty . 

Another possible set up is that of a reflecting wall. 
In this scenario, the grid initially contains a supersonic 
flow, the upstream boundary is set to free-inflow, and the 
downstream boundary is set to a pure reflector (in practice 
this is usually achieved with ghost cells which are exact 
copies of the cells close to the boundary but with their ve- 
locity vector reversed). As the simulation evolves a shock 
is launched from the reflecting boundary and moves up- 
stream into the oncoming flow. Shocked gas which subse- 
quently cools piles up against the boundary (some earlier 
simulations in the literature do not use a "pure" reflecting 
boundary, e.g., Imamura et al.^984). 

Since the compression ratio is relatively small for low 
Mach number shocks (see Table ^1 , simulations with re- 
flecting conditions at the downstream boundary required 
a much larger grid than the simulations with alternative 
set ups. This is because in the latter cases the cold gas 
flows off the grid and the calculation is performed in a 
reference frame where the average position of the shock is 
stationary. In contrast, simulations where the downstream 
boundary is reflecting are performed in a reference frame 
where this boundary is stationary. The width of the CDL 
continuously grows during such calculations as cold gas is 
not allowed to flow off the grid. The growth of the CDL 
pushes the shock upstream into the oncoming flow. The 
normalized speed at which the CDL grows is I/7M 2 , and 
this must be accounted for when calculating the velocity 
of the pre-shock gas in order that a shock with the desired 
Mach number is obtained. 

In Fig. [7| we compare time-space diagrams of simula- 
tions with M = 3 and a = —1.5 and the 3 different set 
ups noted above. In the top panel the overstability is self- 
excited by the small numerical errors which exist from the 
initialization. Linear growth is observed for about the first 
50% of the run, after which the amplitude of the oscilla- 
tions saturate and the system enters the nonlinear regime. 
The elapsed time before the system saturates depends on 
the numerical resolution - finer grids have smaller initial- 
ization errors, so the shock is subject to weaker pertur- 
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bations at startup, and the timescale until saturation is 
longer. 

The middle panel of Fig.[7|repeats the impulsively gen- 
erated simulation shown in the third panel of Fig. [SJ In 
the bottom panel of Fig. we show the evolution of the 
overstability for flow running into a reflecting downstream 
boundary. An overstable shock is immediately generated, 
and the thickness of the CDL increases rapidly with time 
due to the low compression. Unlike in the previous setups, 
shock waves in the CDL are reflected off the downstream 
boundary, and have significant influence on the position 
of the shock early in the simulation. However, the oscil- 
lations become very regular as time progresses, although 
only the first 13 oscillations are shown here. In the simu- 
lations in the upper and middle panels of Fig. [7] the large 
distance to the downstream boundary means that there 
is no mechanism for returning waves transmitted into the 
CDL back into the cooling zone of the shock. 

For ease of comparison, Fig. |H1 shows the time evolu- 
tion of the shock position for the simulations shown in 
Fig. but where we have accounted for the build-up of 
the CDL in the reflecting simulation by subtracting the 
time-averaged growth rate in its width from the position 
of the shock. Broad agreement is seen in the amplitude and 
frequency of the oscillations once the simulations begin to 
settle down, particularly between the upper and middle 
plots. The time required for the amplitude of the oscilla- 
tions to saturate in the top plot is a function of the grid 
resolution, since start-up errors are reduced as the reso- 
lution increases. As already noted, the simulation shown 
in the lower plot becomes very regular as time progresses. 
A power spectrum analysis reveals that the fundamental 
frequency is the same in all 3 simulations to within 1%. 

3.3. The effect of the CDL 

The simulations in the previous section lead us to conclude 
that for low Mach number shocks the particular form of 
the boundary conditions makes little difference to the sta- 
bility. It is interesting to note that the continued growth 
of the CDL does not affect the damping of the oscillations. 

The dynamics of the CDL was extensively investigated 
by Walder & Folini l|1996|) . who concluded that the CDL 
can strongly influence the stability of the cooling zone 
(see also Strickland & Blondin I1995[l . The oscillation of 
the CDL introduces time-dependent conditions at the rear 
boundary of the cooling layer, which may then affect the 
evolution of the radiative shock. Sound waves and shocks 
within the CDL (see Fig.[7|) are created by the pressure on 
the upstream side of the CDL varying as the cooling region 
oscillates through its cycle of overstability. The boundary 
between the CDL and the cooling region oscillates in anti- 
phase with the shock because as the cooling region begins 
to lose pressure support, the CDL finds itself overpres- 
surized compared to the upstream flow (and vice- versa). 
These anti-phase oscillations tend to damp the overstabil- 
ity of the cooling region, since movement of the CDL into 



the cooling region helps to maintain the pressure in the 
cooling region, and vice-versa. 

In those cases where the cooling layer disappears at 
the end of each oscillation cycle, the behaviour of the sys- 
tem should approximately be given by the solution of a 
Riemann problem, with the "left" and "right" states the 
CDL and the supersonic upstream flow, respectively. In 
Table0]we list the flow quantities for the "left" and "right" 
states of the impulsive shock generation simulations shown 
in Fig. El and the Riemann solution for the wave speeds 
moving into these states. In each of the Riemann solutions 
noted in Table^J the left and right waves are both shocks, 
and their velocity and that of the contact discontinuity are 
in good agreement with those apparent in the numerical 
simulations in Fig.El The velocity of the resolved state in 
the Riemann solution, w c d, which corresponds to the ve- 
locity of the boundary between the CDL and the cooling 
layer, is negative for both the M = 5 and M = 10 flows, 
and means that this boundary initially moves downstream 
as the cooling layer grows, as shown in Fig.El We also see 
that v c d becomes less negative with increasing M, again 
in agreement with Fig. [5] Finally, the velocity of both the 
left and right waves, v\ and v r , increase as M increases, 
the former being consistent with the shock oscillation am- 
plitude increasing with M. 

While the initial behaviour of the boundary between 
the CDL and the cooling region can be understood in 
terms of the Riemann problem, its later evolution is de- 
pendent on the behaviour of the cooling region, and vice- 
versa. In Fig. El we see that as the Mach number increases 
the oscillation amplitude of the CDL boundary decreases 
and the amplitude of the overstability increases (when 
specified in units of L c ). Increased damping of the oversta- 
bility occurs when the oscillation amplitude of the CDL 
boundary increases, and such feedback between the CDL 
and the cooling gas is clearly dependent on the Mach num- 
ber of the shock. In the general case (i.e. when the cooling 
region does not necessarily completely disappear between 
oscillations), the degree of feedback between the CDL and 
the cooling gas can be quantified by considering how much 
wave energy from the oscillating post-shock flow is trans- 
mitted into the CDL, and how much is reflected. We as- 
sume that to first order we can represent the boundary 
between the cooling post-shock gas and the CDL as a dis- 
continuous boundary between two states with densities p s 
(the post-shock density) and p c (the density of the CDL). 
The reflection coefficient is then (Landau & Lifshitz Tl 984|l 

R= ( P^~P^ )\ (5) 

\PcCc + PsC s J 

where c c and c s are the sound speeds in the CDL and 
immediately post-shock, respectively. Evaluating Eq. [5] 
reveals that R = 0.074, 0.11, 0.20, 0.36, and 0.59 for 
M = 1.4, 2, 3, 5, and 10. Thus, as the Mach number 
increases, an increasing percentage of the incident wave 
energy is reflected back into the cooling region. This is 
consistent with a decreased oscillation amplitude with in- 
creasing Mach number of the boundary between the CDL 
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Table 4. The initial behaviour of the M = 5 and M = 10 simulations shown in Fig. [5] can be determined from the 
solution to a Riemann problem, which yields the resulting wave velocities in the system (v\, w c( j, and w r being the 
velocity into the left state, the velocity of the contact discontinuity, and the velocity into the right state). The left 
state is defined as the CDL, while the right state is defined as the pre-shock flow (see Fig. EJ . While the values quoted 
are for the initial conditions, they should be a good representation of the conditions which exist at the exact moment 
that the cooling layer disappears when a « a cr . 





Left state 




Right state 


Wave velocities 


M 


P V 


V 


p p V 


Vl Vcd v r 


5 


41.67 1.0 - 


0.024 


1.0 0.024 -1.0 


-0.243 -0.052 0.295 


10 


166.67 1.0 - 


-0.006 


1.0 0.006 -1.0 


-0.117 -0.022 0.312 



and cooling zone (see Fig. EJ), since less energy is trans- 
mitted into it. It is clear, therefore, that the CDL is more 
able to damp the shock overstability when the Mach num- 
ber is small, since it is then more efficient at absorbing 
the energy of the oscillations. This efficiency is indepen- 
dent of the thickness of the CDL, and explains why we do 
not see more efficient damping, indicated by a decrease in 
the oscillation amplitude of the cooling region, when the 
thickness of the CDL grows (see Figs. and EJ. At high 
Mach numbers the fraction of transmitted wave energy ap- 
proaches zero, the velocity of the CDL boundary and the 
amplitude of its oscillations become very small, and feed- 
back between the CDL and the cooling region vanishes. 
The Mach number dependent transmission efficiency of 
sound wave energy is undoubtedly the reason why we ob- 
serve a strong Mach number dependence of a cr and U{ at 
low Mach numbers, but not at high Mach numbers (where 
R reaches an asymptotic value of 1.0). 

Our results are consistent with those presented by 
Walder & Folini H199fi[) . who noted from their studies 
of shocks with higher Mach numbers that feedback be- 
tween the CDL and the cooling region is relatively mi- 
nor if the pre-shock gas is smooth (in agreement with our 
M = 10 simulations). They also found that the feedback 
is stronger if the shock encounters a disturbance in the 
pre-shock medium. In such situations the CDL may os- 
cillate with significantly greater amplitude, which in turn 
may cause the overstability to increase or decrease in am- 
plitude and/or frequency, and even switch on and switch 
off (see Figs. 15 and 16 in Walder & Folini HMfy . While 
we have not explored the effect of such upstream distur- 
bances in this work, we anticipate that these changes in 
behaviour will apply equally to low Mach number shocks. 

An interesting point in our simulations is that energy 
transmitted into the CDL can leave it in two different 
ways. Waves in the CDL will be damped through radia- 
tive losses as compressions increase the temperature above 
To. In addition, energy may be lost from the system if the 
downstream boundary on the hydrodynamical grid has 
outflow conditions (though no such losses occur in the sim- 
ulations where we impose a reflecting wall at the back of 
the CDL). In the simulations noted in Sec. 13. II the down- 
stream boundary is far enough away that the transmitted 
waves never reach it. 



We also note that as the CDL is thick behind low Mach 
number shocks, such shocks should be less susceptible to 
other types of instability than their higher Mach number 
counterparts, such as the non-linear thin shell instability 
(Vishniac HHMI Blondin & Marks HMty . and non-radial 
instabilities. 

3.4. Cooling below the pre-shock temperature 

In all of the simulations in the previous sections, the 
shocked gas is prevented from cooling below the temper- 
ature of the pre-shock gas. Since there is no fundamental 
reason why this could not occur, we have examined the 
nature of the overstability when x = T c /Tq < 1. 

We find that the value of x has a dramatic effect on 
the properties of low Mach number shocks. Fig. HOIshows 
the time evolution of the position of the shock front as a 
function of x when M — 1.4 and a = —2.5. It is obvi- 
ous that decreasing x causes the oscillation amplitude to 
increase and the frequency of the fundamental mode to 
decrease, which is the same trend obtained when M in- 
creases with x = 1 and ol fixed. Clearly the value of x has 
a strong influence on the behaviour of the shock, though 
we see convergence at low values of x (the behaviour of 
a simulation with M = 1.4, a = —2.5, and x = 10 -3 , is 
not too different from that with x = 10~ 2 ). At high Mach 
numbers (e.g., M = 40), varying x has a negligible effect. 
Fig. El shows that the critical value of a for stability is 
also dependent on x- In the top panel the overstability is 
damped, but this is not the case in the lower panel. In 
Table we note values of a cr as a function of both M 
and x- There are large changes in a cr with x at low Mach 
numbers, with the result that it is much easier for a low 
Mach number shock to be overstable when x < 1- We do 
not obtain a cr = 0.4 at lower Mach numbers when x is 
reduced. 

It is clear from Figs. 1101 andll Hand Table that reduc- 
ing x below unity makes the shock behave like a higher 
Mach number shock with x = 1- This is expected since if 
X < 1, the density of the CDL increases above the value 
given by 7M 2 (the cooling time and length of the shocked 
gas also increase, though by amounts which are dependent 
on a). A shock with M — 1.4 and x — 0-1 has p c /pa = 42, 
which is almost identical to a shock with M = 5 and x = 1 • 
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Table 5. Critical value of a for damping of the funda- 
mental mode as a function of the Mach number and Xi 
the ratio of the temperature of the CDL to the pre-shock 
temperature. This table is an extension of Table |3] where 
a cr is listed for x = 1- 



M 






X 






1.0 


0.5 


0.1 


0.01 


1.40 


-1.8 


-0.6 


-0.2 


-0.1 


2.00 


-1.2 


-0.6 


-0.2 


-0.1 


3.00 


-0.7 


-0.5 


-0.2 


-0.1 


5.00 


-0.4 


-0.3 


-0.1 


-0.1 


10.0 


-0.1 


-0.1 


0.1 


0.3 


20.0 


0.2 


0.2 


0.2 


0.3 


40.0 


0.3 


0.3 


0.3 


0.3 


100 


0.4 


0.4 


0.4 


0.4 



However, the value of the reflection coefficient, R — 0.57, 
is almost equivalent to that from a shock with M = 10 
and x — 1- It is n °t surprising, therefore, that the value of 
a CI for a shock with M = 1.4 and x — 0.1 is in good agree- 
ment with the value obtained for a shock with M = 10 
and x — 1< Thus low Mach number shocks with x < 1 are 
more susceptible to overstability than a straightforward 
comparison of p c / Po would suggest. An evaluation of R in 
this case allows a much more accurate estimate of a CT . 

Values of R as a function of M and x are shown in 
Table In Fig. ^] we display the value of ev cr as a func- 
tion of R for 4 different values of x- It is clear from the 
relatively tight relation between the curves that the value 
of R can be used to obtain a reasonable estimate of ev cr . 
The larger scatter in a cr for a given R which exists in 
the top right of the plot indicates that our simple method 
of evaluating R is a somewhat poorer description of the 
actual dynamics in this region of parameter space. 

We also find that not all of the shock properties scale 
in such a straightforward manner. For example, the oscil- 
lation frequency of the fundamental mode when M = 1.4, 
a = —2.5, and x — 0.1 is uif = 0.31, which is some 
way from the value Wf = 0.14 obtained when M = 10, 
a = —2.5, and x — 1 ( m comparison, when M — 1.4, 
a = —2.5, and % = 1, we find u)t — 0.67). Thus the actual 
Mach number of the shock is still important. 

4. Discussion 

Radiative shocks occur ubiquitously in astronomy, e.g., in 
SNRs, stellar wind bubbles, and jets from young-stellar- 
objects. They exist over a wide range of spatial scales, 
from accreting white dwarfs to galactic scale-accretion. 
Evidence consistent with the presence of an unsteady 
flow, such as produced by the radiative overstability, has 
been found in observations of the Vela SNR (Raymond et 
al. H991(l - unusually broad line widths (of Sin and Mgn) 
and differences between the shock ram pressure and the 
postshock thermal pressure. The use of steady shock mod- 
els in such circumstances means that conclusions drawn 



Table 6. Reflection coefficient, R, as a function of the 
Mach number and Xi the ratio of the temperature of the 
CDL to the pre-shock temperature. 



M 






X 






1.0 


0.5 


0.1 


0.01 


1.40 


0.075 


0.239 


0.567 


0.839 


2.00 


0.115 


0.261 


0.575 


0.842 


3.00 


0.199 


0.342 


0.632 


0.866 


5.00 


0.356 


0.490 


0.731 


0.906 


10.0 


0.588 


0.688 


0.847 


0.949 


20.0 


0.765 


0.828 


0.919 


0.974 


40.0 


0.875 


0.910 


0.959 


0.987 


100 


0.948 


0.963 


0.983 


0.995 



from the observed line ratios will be incorrect (fnnes et 
al. H987al [l987bl Gaetz et al. lT3S£|) . 

A key point concerning SNRs is that the deceleration 
of the blast wave limits the number of oscillations that 
may occur, since there is only a finite time during which 
the shell is unstable (see, e.g., Mac Low & Norman 1993 
Blondin et al. 1998). As noted in Sec. □ previous works 
have assumed that the temperature of the ISM corre- 
sponds to the WIM phase (i.e. T 1SM « 10 4 K; see McKee & 
Ostriker 1977). The consequence of this is that low Mach 
number shocks have post-shock temperatures which are 
on the steeply rising part of the cooling function, and will 
thus be stable. However, the ISM is known to consist of 
a variety of phases, with its volume likely dominated by 
gas with a temperature of ~ 10 5 — 10 6 K (the HIM phase 
in the nomenclature of McKee & Ostriker ll977fl . This hot 
phase is continuously regenerated through the combined 
actions of supernova explosions and the winds of massive 
stars. SNRs interacting with this phase of the ISM have 
much reduced Mach numbers. 

The evolution of SNRs can be broken up into a se- 
quence of stages. From an initial ejecta-dominated stage, 
they may evolve through the Scdov- Taylor (ST), pressure- 
driven snowplough (PDS), and momentum-driven snow- 
plough (MDS) stages. Once the SNR is in the PDS stage, 
its forward shock may be susceptible to the radiative over- 
stability. The timescale for the transition between the ST 
and PDS stages is roughly (Blondin et al. H998fl : 

i rad w 2.9 x 10 4 i4 /17 n 9/17 yr, (6) 

where E^i is the kinetic energy of the original explosion 
in units of 10 51 ergs, and uq is the number density of 
the ISM. The velocity of a ST blast wave at this age is 
v s w 260£ , 5{ 17 ?1q/ 17 kms -1 . The sound speed in the hot 
ISM is c w 1501^,6 kins" 1 , where Ti SM ,6 is the ISM 
temperature in units of 10 6 K. 

For the blast wave to be susceptible to overstability 
when the SNR is expanding at relatively low Mach num- 
bers requires that the post-shock gas is located on a part 
of the cooling curve where a < a CI . The simulations in 
Sec. 01 show that a < —0.7, —1.2, and —1.8 is needed in 
order to have overstable shocks at Mach numbers of 3, 
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Table 7. Approximate properties of SNR blast waves at 
the transition into the pressure-driven snowplough stage 
for different values of Tjsm and n , chosen to demon- 
strate the conditions under which SNRs in the pressure- 
driven snowplough stage have low Mach number forward 
shocks which are susceptible to the radiative overstabil- 
ity. The age of the SNR at the onset of shell formation, 
tn, the shock velocity, v s , and the Mach number, M, at 
this time, are given. The values illustrate that overstable 
low Mach number radiative shocks in SNRs will exist. For 
Xism = 10 6 K and n = 3.5 x 10~ 3 cm~ 3 the blast wave 
becomes subsonic before transition to the pressure-driven 
snowplough stage. 



2ism 


no 


tti 


v B 


Ts/Tism 


M 


(K) 


(cm" 3 ) 


(10 5 yr) 


(kms^ 1 ) 






2 x 10 s 


10 -3 


11.2 


115 


1.7 


1.7 




1(T 2 


3.3 


151 


2.4 


2.3 




lcr 1 


1.0 


198 


3.6 


3.0 




l 


0.3 


260 


5.6 


3.9 


10 6 


10 -3 


11.2 


115 




0.8 




1(T 2 


3.3 


151 


1.0 


1.0 




KT 1 


1.0 


198 


1.5 


1.5 




1 


0.3 


260 


1.7 


1.7 



2, and 1.4, respectively. Such values of a cr are fairly de- 
manding in the sense that a is less than a cr only in very 
limited temperature ranges. Nonetheless, SNRs evolving 
into an ISM with Ti SM ~ 10 5 K or Ti SM > 10 6 K may 
satisfy these conditions (see Fig.|5J). However, the work in 
Sec. l3.4l showed that the criteria for overstability becomes 
far less stringent if the post-shock gas cools to a lower tem- 
perature than the pre-shock gas (a cr « —0.1 when x = 0.1 
for M > 1.4), in which case SNR shocks are readily over- 
stable. Table shows values for t ra d, v s , T s and M for two 
values of Ttsm and for different values of no, assuming 
-E51 = 1 and 7 = 5/3. Since itr is an approximation to the 
time of shell formation the numbers in TableQare only il- 
lustrative. Nevertheless, Table0shows that SNRs evolving 
in high pressure environments may have low Mach number 
shocks which are susceptible to the radiative overstability, 
as determined from the slope of the cooling curve between 
Tism and T s and by the Mach number dependence of a cr 
(see Table EJ). Our work in this paper is therefore relevant 
to such extreme settings as starburst galaxies, where the 
range of density and temperature values in Table \7\ are 
comparable to observed ranges. 

It has recently become apparent that radiative shocks 
in dense molecular clouds can also be unstable (Smith & 
RosenEHnni- While the simulat ions presented in this work 
are of relatively high Mach number shocks, slower shocks 
with Mach numbers of about 5 and post-shock tempera- 
tures of about 10 3 K, should also be overstable (see Fig. 2a 
in Smith & Rosen 2003), though as usual there is the 
caveat that magnetic fields tend to stabilize the oversta- 
bility (Toth & DraineHnSI- 



5. Summary 

In this paper we have examined the stability of low Mach 
number radiative shocks. We find that such shocks may be 
overstable for sufficiently low values of a, and that there 
are temperature regions in representative cooling curves 
where the required values of a are obtained. We have de- 
termined the critical value of a at the boundary between 
stability and overstability as a function of Mach number, 
finding that a cr increases with M. The strong shock limit 
of a cr « 0.4 for damping of the fundamental mode is not 
reached until M « 100. The frequency and amplitude of 
the fundamental mode of oscillation are strongly depen- 
dent on the Mach number of low Mach number shocks. 
Feedback between the cooling region and the CDL is also 
a function of the Mach number, and may be quantified in 
terms of the reflection coefficient of sound waves in the 
general case. In those cases where the cooling layer com- 
pletely disappears at the end of each oscillation cycle, the 
velocity of the shocks driven upstream into the pre-shock 
flow and downstream into the CDL, and the velocity of 
the boundary between the CDL and the cooling layer, are 
in good agreement with those determined from a solution 
to the Riemann problem. 

An important finding is that the stability properties of 
low Mach number shocks are dramatically altered if the 
shocked gas is able to cool below the temperature of the 
pre-shock gas. In such circumstances, low Mach number 
shocks behave in many ways as if they had higher Mach 
numbers. An increase in cv cr is the most fundamental of 
these changes, meaning that such shocks are more likely to 
be overstable. An investigation into the effects of differing 
conditions for the initial set up and the grid boundaries 
reveals that these have very little influence on the stability 
criteria and oscillation frequencies. This is due to the fact 
that in low Mach number shocks the CDL acts much more 
like a cushion than a reflector, absorbing a substantial 
amount of the incident sound wave energy. Growth in the 
thickness of the CDL does not enhance the damping of 
the oscillations. 

Our work is relevant to the evolution of SNRs inter- 
acting with the hot phase of the ISM (e.g., in starbursts), 
and to shocks in dense molecular clouds. 
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Fig. 6. Power spectra (mean square amplitude) for the 
simulations shown in Fig. where a = —1.5. M varies 
from 1.4, 2, 3, 5 from the top to bottom panels respec- 
tively. 
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Fig. 7. Time-space diagrams of the density evolution of 
a 1-dimensional radiative shock with M = 3 and cooling 
exponent a — —1.5. In the top panel the flow is initial- 
ized to the steady state solution. In the middle panel the 
shock is impulsively generated. In the bottom panel the 
downstream boundary condition is reflecting. 
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Fig. 8. The shock position as a function of time for simu- 
lations with M = 3 and a = —1.5. The set up with an iso- 
lated steady-state shock as the initial condition is shown 
at the top, the set up with impulsive shock generation 
is shown in the middle, and the set up with a reflecting 
downstream boundary (i.e. "flow into a wall" ) is shown at 
the bottom. In the latter case the growth in width of the 
cold dense layer has been removed to aid comparison with 
the other scenarios, and each datasct is offset with respect 
to the others. 
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Fig. 9. Time-space diagrams of the density evolution of a 
1-dimensional radiative shock with cooling exponent a = 
-1.5 and M = 5 (top) and M = 10 (bottom). 
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Fig. 10. The time-evolution of the shock position as a 
function of \ when M = 1.4 and a = —2.5. 
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Fig. 11. The shock position as a function of time when 
M = 1.4 and a = —1.5. The ratio of the temperature of 
the CDL to the pre-shock temperature, x, is noted in each 
panel. 
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Fig. 12. The value of a cr , the critical value for damping of 
the overstability, as a function of the reflection coefficient, 
R, for four different values of \- The lines are only meant 
to guide the eye. 



